Model ab initio study of charge carrier solvation and large polaron formation on 

conjugated carbon chains 
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Using long CjvEb conjugated carbon chains with the polyynic structure as prototypical examples 
of one-dimensional (ID) semiconductors, we discuss self-localization of excess charge carriers into 
ID large polarons in the presence of the interaction with a surrounding polar solvent. The solvation 
mechanism of self-trapping is different from the polaron formation due to coupling with bond-length 
modulations of the underlying atomic lattice well-known in conjugated polymers. Model ab initio 
computations employing the hybrid B3LYP density functional in conjunction with the polarizable 
continuum model are carried out demonstrating the formation of both electron- and hole-polarons. 
Polarons can emerge entirely due to solvation but even larger degrees of charge localization occur 
when accompanied by atomic displacements. 

PACS numbers: 31.15.A-, 31.15at, 71. 15. Mb, 71.38.Fp, 81.07.Nb 



I. INTRODUCTION 



One-dimensional (ID) semiconductor (SC) struc- 
tures such as 7T— conjugated polymers, nanotubes and 
nanowires are interesting nanoscopic objects that can 
be exploited in various areas, including (opto)electronics, 
energy harvesting and sensors. The nature and proper- 
ties of excess charge carriers on these structures are fun- 
damental for many applications. As covalent bonds pro- 
vide for wide electronic bands and effective band masses 
much smaller than the free electron mass, intrinsically ex- 
cess charge carriers in these systems have a high propen- 
sity for delocalization and a potential for high mobili- 
ties. If not for extrinsic defects, inherent limitations on 
the mobility arise due to interactions of charge carri- 
ers with other subsystems - most commonly with dis- 
placements of the underlying atomic lattice of the ID 
SC structure. This electron-phonon coupling has been a 
subject of numerous studies for various specific ID sys- 
tems. It is well-known that, as a result of the inter- 
action with lattice phonons, excess charge carriers may 
undergo self-localization into polaronic states. Such po- 
larons have been extensively discussed in the context 
of conjugated polymers (CPs), where they are believed 
to be accompanied by localized bond-length modulation 
patterns and new features in the optical absorption (see, 
e.g., Refs. [3, for reviews and original references). 

A different and much less explored implementation 
of the strong polaronic effect can take place when a 
ID semiconductor is immersed in a 3D polar medium, 
the situation of particular relevance to applications and 
processes involving fundamental redox reactions in po- 
lar solvents. In this case of excess charge carrier solva- 
tion, the long-range Coulomb interaction, as has been 
recently discussed within simplified theoretical models 
[3> 0) B B 0]' can result in the formation of ID adia- 
batic large-radius polarons, where a localized electronic 
state on a SC structure is surrounded by a self-consistent 
pattern of the sluggish (orientational) polarization of the 



solvent. The Coulomb mechanism of the polaron forma- 
tion here is analogous to the well-known 3D polarons in 
polar SCs [EM E3, [H| and solvated electrons in polar 
liquids [IiI.1pX|13||. Distinct from those cases are the 
confinement of the electron motion in a ID nanostruc- 
ture and its structural separation from the 3D polariz- 
able medium. The energetic significance of the solvation 
has been emphasized by finding 0, H[ that the binding 
energy of the resulting ID polarons could reach a sub- 
stantial fraction, roughly one third, of the binding en- 
ergy of Wannier-Mott excitons, the well-known primary 
photoexcitations in many ID SCs. This may lead to en- 
hanced charge separation. On the other hand, the mobil- 
ity of solvated charge carriers is drastically reduced due 
to the dissipative drag of the medium 0, [y, [l4[ • New op- 
tical absorption signatures are also expected as a result 
of solvation @,[H- 

Given the generic character of the solvation mecha- 
nism of charge-carrier self-localization and its potential 
implications for carrier mobilities and redox processes on 
ID SC nanostructures, this effect appears interesting and 
important enough to warrant studies at various compu- 
tational levels. As a step in that direction, we have 
attempted a first - to our knowledge - model ab ini- 
tio study reported in a recent communication [lfj that 
would demonstrate the formation of solvation-induced 
ID large-radius polarons. In this paper we elaborate on 
that brief report discussing features of solvation-induced 
self-localization of excess charges on finite linear conju- 
gated carbon chains (polyyne oligomers) as de- 
rived within density-functional theory (DFT) compu- 
tations in conjunction with the polarizable continuum 
model (PCM) 0, [HI, 01 . 

As in numerous studies of electron-lattice polarons 
in CPs, ab initio studies are expected to be helpful 
in elucidating the role played by valence electrons and 
many-electron interactions in accommodation of an ex- 
tra charge carrier. It should be noted that a reorgani- 
zation of valence band electrons can be important even 
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in non-interacting electron models - a nice illustration 
was given in Ref. [l9| that analyzed how polarons of 
a two-band Peierls dielectric model evolve into single- 
particle Holstein (2(| polarons in the limit of the "frozen 
valence band" approximation. With realistic Coulomb 
interactions in place, the role of valence electrons in 
the formation of the relevant self-consistent potentials 
would only increase. Interestingly enough, despite a rel- 
atively long history of first-principles studies of electron- 
phonon polarons in CPs (see, e.g., multiple references 
in Refs. [2l[ I22!]). certain questions have been raised re- 
cently regarding applicability of different ab initio frame- 
works to describe the formation of those polarons. That 
pertains to the failure of the local-density-approximation 
and generalized-gradient-approximation DFT schemes to 
detect self-localized charge density distributions in var- 
ious charged oligomers, while Hartee-Fock, parameter- 
ized semi-empirical and possibly hybrid-functional DFT 
methods have been reported to lead to c harg e localization 
in the middle of oligomers (e.g., Refs. [HlUE 53, HH)- 
Quite different level-of-theory-dependent magnitudes of 
the effective electron-phonon coupling have also been 
found even in neutral polymeric systems (see, e.g., a com- 
parative discussion in Refs. [H,[23|). A caution therefore 
has to be exercised in interpretation of at least the quan- 
titative aspects of various ab initio results. 

As first-principles calculations are applied to the ef- 
fect of solvation on charge-carrier self-localization within 
the PCM framework, it would also be desirable to have 
a comprehensive comparison of different computational 
approaches. Such a study is currently underway and we 
are planning to report its results elsewhere. Reassur- 
ingly, the effect of solvation on excess charge localiza- 
tion appears qualitatively robust in those computations. 
While providing an illustrative example of comparison 
with Hartree-Fock computations, in this paper we rather 
focus on a more modest goal to discuss trends that are 
derived within a single computational scheme. Specifi- 
cally, we employ DFT calculations with the B3LYP hy- 
brid exchange-correlation density functional [28|. This 
functional is known to include both local and non-local 
effects and has been commonly used in recent studies 
EE HE [H [H of polyyne and its oligomers. 

By exploring model systems with prescribed carbon 
atom positions, we will be able to illustrate excess charge 
carrier (an electron or a hole) localization and the re- 
sulting formation of large ID polarons entirely due to 
the polarization of the surrounding solvent. Examples of 
self-consistent - that is, practically independent of the 
chain length - polaron structures will thereby be also 
given. Within the confines of the computational scheme 
used, "realistic" geometries of carbon atom arrangements 
are those that optimize the total energy of the system. 
We have performed such geometrical optimizations for 
polyynic chains both in vacuum and in solvent. In the 
range of chain lengths studied (up to TV = 100 of carbon 
atoms) with the B3LYP hybrid functional, we have not 
been able to detect a clear trend towards excess charge 



self-localization due to the electron-phonon interaction 
alone - as the chain length increased, the bond-length 
modulation pattern would only become flatter in the mid- 
dle of the chain. On the contrary, in the presence of the 
Coulomb interaction with the solvent, clearly localized 
patterns of bond-length alternation have been found to 
act as to boost the degree of excess charge localization, 
possibly reflecting synergistic effects from both mecha- 
nisms in this case. 



II. SYSTEMS, COMPUTATIONS, AND DATA 
ANALYSIS 

Since the computational demand increases substan- 
tially in DFT-PCM calculations, we have selected struc- 
turally simple ID semiconducting systems of hydrogen- 
terminated polyynic linear carbon chains CjvH2 for our 
demonstration. It should be stressed however that, while 
playing the role of a prototypical example in our study, 
polyynic chains continue to be the subject of much atten- 
tion in their own right ( [H, HE HH, H2] and multiple refer- 
ences therein) . Interestingly, they were predicted [H, H3] 
to possess a rich family of electron-lattice self-localized 
excitations because of the extra degeneracy of molecular 
orbital (MO) levels. 

The subject of our attention in this paper are long 
chains CatH2 with even number of carbons TV (TV between 
20 and 100). It is well established now [H H3, HI, HI, [H| 
that even-TV long (TV > 10) neutral chains have their 
ground state with the polyynic structure, that is, they 
feature a (nearly uniform) alternating pattern of triple 
C = C and single C — C bonds leading to a gap in the 
electronic spectrum. The ground state of odd- TV chains 
would exhibit a kink-like "defect" in the bond-alternation 
pattern with associated mid-gap states [3(1 HH, which 
we do not discuss here. Having neutral chains CatH^ as 
a benchmark, we have studied both negatively charged 
CatH^ (an extra electron) and positively charged CatH^ 
(an extra hole) chains. As is clearly illustrated below, 
their self-localization behavior has been found very sim- 
ilar, with the excess charge measured correspondingly in 
units of (— e) or (e), where e is the magnitude of the 
fundamental charge. To save space, some of the results, 
therefore, would be shown only for the negatively charged 
systems. 

For the ab initio engine in this study we employ den- 
sity functional theory (DFT) within the Gaussian 03 [l6[] 
set of programs. More specifically, the majority of the re- 
sults presented in this paper have been derived with the 
B3LYP hybrid exchange-correlation functional which in- 
cludes the Hartree-Fock exchange along with DFT corre- 
lation [28| . In order to investigate the effects of solvation 
in a polar medium, Gaussian 03 [l6| offers its implemen- 
tation of the Polarizable Continuum Model (PCM) de- 
scribed in original publications [13, [H, [H, S3] ■ A rich 
6-311+- |-G (d,p) all electron basis set has been employed 
throughout. It should be noted that the results of our "in 
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FIG. 1: Spatial charge distributions on N — 100 polyynic chains with RG1 geometry derived by using atomic charges from 
Mulliken (panel columns indexed 1 and 3) and Lowdin (panel columns 2 and 4) procedures. Panel row (a) is for the totally 
neutral chains, row (b) for the negatively charged chains with one added electron; row (c) depicts the distribution of excess 
charge defined as a difference of the data in rows (b) and (a). Two left columns (1 and 2) depict charges per atom, two right 
columns (3 and 4) yield charges per unit cell. In each of the panels, solid lines show the results for chains in solvent, whereas 
dashed lines (when distinguishable) for chains in vacuum. It is evident that cell-centric excess charge distributions in panels 
(c3) and (c4) are nearly identical despite substantial differences between Mulliken and Lowdin charges observed in other panel 
pairs. 



vacuum" calculations have been verified to compare well 
against previously published ab initio data on neutral [27} 
and charged polyynic systems both in terms of en- 
ergetics and optimized bond-lengths. For all DFT-PCM 
("in solvent") calculations, water has been chosen as a 
polar solvent with its default parameters in Gaussian 03. 
We remind the reader that the focus of this paper is to 
compare results obtained in vacuum and in solvent envi- 
ronments within the same computational scheme rather 
than compare various computational approaches. 

In order to be able to clearly discern the effects due 
to the solvation and provide conceptual proofs and illus- 
trations (see Sec. Mil) , we will be exploring both systems 
with the full geometric optimization of the underlying 
atomic lattice and systems with model rigid geometries, 
that is, with prescribed atomic positions. One of such 
rigid geometries, referred to as RG1, features the follow- 
ing fixed lengths: 1.235 A for the triple bonds, 1.326 A 
for the single, and 1.063 A for the C — H bonds. These 



lengths have been retrieved from our complete optimiza- 
tion of the neutral N — 80 chain in vacuum and are 
in fact relatively close to optimal values (see Sec. IIVI) . 
Other rigid geometries have also been explored where all 
single C — C bonds are artificially stretched in successive 
increments of 0.1 A up to the final length of 1.726 A, 
the latter geometry is denoted as RG2. This artificial 
stretching is used to study the effect of narrower elec- 
tronic bandwidths (larger band effective masses) on the 
charge localization. 

A very relevant quantity for our discussion is the equi- 
librium spatial charge distribution over the ID SC, which 
is naturally related to atomic charges in outputs of ab 
initio calculations. Importantly, these charges are cal- 
culated from many-electron wave functions and thereby 
reflect responses due to all electrons in the system. Two 
procedures, due to Mulliken and Lowdin, are widely used 
for charge population analysis [13] • It is well known [ItJ 
that calculations of atomic charges depend on the ba- 



sis set transformations and sometimes lead to artifacts 
and spurious results. Such artifacts are in fact appar- 
ent in our illustration in Fig. [1] that compares various 
results we obtained for N = 100 chains. Differences be- 
tween Mulliken and Lowdin charges are especially clearly 
seen for raw charge populations separately on neutral 
(row (a)) and charged (row (b)) chains, whether they 
are atomic- or unit-cell-centric. Comparison of panels 
(cl) and (c2), however, indicates that certain cancela- 
tion effects take place and both methods lead to less 
dissimilar results when one is interested in the spatial 
distribution of the excess charge derived as the differ- 
ence of charge densities on charged and neutral chains. 
Still those atomic-centric results exhibit spurious oscilla- 
tions which may be artifacts related to the presence of 
bond-oriented charge-density waves in our system. This 
suggests "averaging" of the charge over the unit cells con- 
sisting of pairs of neighboring carbons (for convenience, 
the chain end cells are defined here to include the hydro- 
gen atomic charges). Panels (c3) and (c4) show that var- 
ious artifacts indeed practically disappear, and the unit- 
cell-based excess charge distributions calculated by Mul- 
liken and Lowdin procedures become remarkably close. 
We confirmed this conclusion for many other cases of in- 
terest. Accordingly, it is these stable results for excess 
charge per cell that will be displayed in our results below 
(using raw Lowdin charges). 



III. SOLVATION-INDUCED SELF-CONSISTENT 
CHARGE LOCALIZATION IN MODEL SYSTEMS 

In real systems, both atomic displacements and po- 
larization of the surrounding medium are expected to 
take place in order to accommodate an excess charge 
carrier. Moreover, given the nonlinear nature of the po- 
laronic effect, both subsystems may act in a synergistic 
way. In order to better understand the contribution com- 
ing from the solvation, we therefore start from studying 
model systems with rigid geometries, where no localized 
bond-alternation patterns are allowed while the medium 
"adjusts" its state of polarization. 



A. Spatial distribution of the excess charge 

Figures [2] and [3] compare excess charge density distri- 
butions, negative and positive, respectively, calculated 
as described in Sec. UH for chains of varying lengths both 
in vacuum and in solvent. It is evident that panels 
(a) and (b) of each of the figures display quite distinct 
trends. Disregarding the end effects (end cells include 
charges on hydrogens), vacuum results in panels (a) re- 
flect the "particle in a box" behavior: as the chain length 
increases, the excess charge is distributed more and more 
uniformly over the whole chain. In a sharp contrast, for 
chains in the solvent, panels (b), the excess charge ex- 
hibits much more localized distributions around central 
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FIG. 2: Excess charge density spatial distribution on nega- 
tively charged CjvHj chains with N = 20, 40, 60, 80 and 100 
- top to bottom curves, respectively, in their central parts. 
All chains feature the same pattern of fixed bond-lengths cor- 
responding to RG1 geometry described in text. Panel (a) 
displays results for chains in vacuum, panel (b) for chains in 
the solvent environment. 
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FIG. 3: As in Fig. [2] but for positively charged CatHJ chains 
with N = 20, 40, 60 and 80 - top to bottom curves, respec- 
tively, in their central parts. 

parts of the chains, to the extent that no appreciable end 
effects are present. As the chain length increases, those 
distributions show a very clear tendency towards conver- 
gence, albeit not completely achieved within the range of 
lengths studied. 

Instead of further increasing the chain length in al- 
ready demanding computations, we now attempt to de- 
crease the spatial extent of charge localization in model 
systems. To this end, we recall that, in single-particle 
models 0, 0, [E IE > the spatial extent of the polaron 
is determined by the balance of the gain in the poten- 
tial energy of the self-localized carrier due to its interac- 
tion with the medium, and the loss in its kinetic energy 
due to localization. Decreasing the role of the kinetic en- 
ergy would lead to a diminished derealization propensity 
and is expected to shorten the polaron. This should be 
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FIG. 4: Excess charge density spatial distribution on charged 
CjvHjT chains, (a) N=80. The evolution of the charge distri- 
bution upon stretching the single bonds by 0.1 A successive 
increments from 1.326 A to 1.726 A - bottom to top in a 
group of 5 clearly distinct curves shown by different connected 
symbols and denoted "in solvent". The same - but uncon- 
nected - data symbols are used to show the results for chains 
in vacuum, those distributions largely overlap, (b) Compar- 
ison of distributions for chains in the solvent with N = 40, 
60, 80 (different data symbols) and 100 (solid line) for the 
geometry RG2 featuring the longest single bonds. Curves for 
N = 60, 80 and 100 practically coincide. 
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FIG. 5: Excess charge density spatial distributions on 
charged CjvHjJ" chains, (a) N=80. The evolution of the 
charge density distribution upon stretching bonds as in Fig. [4] 
(b) Comparison of distributions for chains in the solvent, 
with N = 40, 60, 80 (different data symbols) for the RG2 
geometry. Curves for TV = 60 and 80 practically coincide. 



achievable with an increased effective band mass of the 
carrier, or with a narrower electron bandwidth. To affect 
this, we chose to artificially stretch all single bonds in 
our rigid model system as described in Sec. [Til thereby 
making C = C "dimers" more and more separated from 
each other. Figures HJa) (extra electron) and [5£a) (ex- 
tra hole) illustrate a dramatic difference in the response 
of the excess charge distributions on N — 80 chains to 



successive single-bond stretchings in vacuum and in sol- 
vent. Whereas the charge density in vacuum shows only 
very little changes in this process and preserves nearly 
uniform distributions over the whole chain, for chains in 
solvent, each successive stretching step indeed results in 
substantially smaller and smaller localization regions. As 
the extent of localization in the longest-stretched geome- 
try, RG2, appears to be well within the chain length, we 
can now compare the corresponding excess charge dis- 
tributions for several chain lengths. Figure 0Jb) makes 
it evident that the distributions for N — 60, 80, and 
100 negatively charged chains in solvent practically coin- 
cide. Likewise, distributions for TV = 60 and 80 positively 
charged chains are practically identical in Fig. [S[b) • This 
confirms that the obtained charge distributions on longer 
chains correspond to a fully converged self-consistent 
pattern of excess charge localization. Importantly, self- 
localization of the excess charge demonstrated here has 
been achieved entirely due to the interaction with the 
surrounding solvent. 

Comparing results for negatively charged and posi- 
tively charged species, Fig. [2] vs Fig. [3] and Fig. [4] vs 
Fig. one, of course, notices some effects of charge- 
conjugation (CC) symmetry breaking. They depend on 
the chain length as affected by both the ends and by the 
inherent band structure. As the inherent 7r-orbital bands 
are nearly CC-symmetric, these effects are evidently rel- 
atively weak, and basically the localization behavior for 
electrons and holes appears very similar. 



B. Molecular orbitals 

We now turn from the total excess charge distribution 
to individual MOs and the resulting electronic structure. 
In the neutral CatH^ chains, 7r-electron energy levels ex- 
hibit a four-fold degeneracy - in addition to the spin de- 
generacy, there is a two-fold degeneracy with respect to 
the spatial orientation of 7r-orbitals. If the chain direc- 
tion is chosen as the z-axis, 7r-orbitals can be said to 
be oriented either along x or along y axes. In charged 
chains, spin- and orbital-orientation-degeneracy are gen- 
erally lifted. Also, in the presence of the polaronic ef- 
fect, one observes both spatially localized and delocalized 
(over the chain length) individual electronic states. The 
overall picture of MO energy levels in neutral and charged 
chains is illustrated in Figs. [5] (CiooHij vs CiooH^) and 
[7] (CsoHj vs CsoH^), while Fig. [5] illustrates the spatial 
behavior of MOs on charged CiooH^" chains in more de- 
tail. 

Looking at the electronic structures of neutral and 
charged chains in Figs. [5] and one can compare the 
differences in responses between chains in vacuum and 
chains in solvent as well as compare the differences caused 
by the geometry of carbon atom arrangements (RG1 vs 
RG2) . Going to the artificial RG2 geometry with longer 
C — C bonds allows us here to decrease the polaron size 
and to illustrate some features that would otherwise be 
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FIG. 6: Molecular orbital energy levels for neutral (0) and 
charged (-1) C100H2 chains with rigid RG1 (a-d) and RG2 (e- 
h) geometries both in vacuum and in the solvent. The arrows 
indicate positions of HOMO levels. 
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FIG. 7: Molecular orbital energy levels for neutral (0) and 
charged (+1) C80H2 chains with rigid RG1 (a-d) and RG2 (e- 
h) geometries both in vacuum and in the solvent. The arrows 
indicate positions of LUMO levels. 



hard to observe within the range of chain lengths stud- 
ied. Probably the first visually noticeable difference be- 
tween RG1 and RG2 electronic structures is a much wider 
band gap in RG2 - as indeed one should expect based on 
a much larger difference in the hopping integrals corre- 
sponding to triple and single bonds. As such, however, 
the size of the band gap is inconsequential for the de- 
gree of charge localization due to solvation. Instead, of 
direct importance here in a much narrower bandwidth 
in RG2, which is easily detected via a higher density of 
states (panels (e) vs (a) and (g) vs (c)). The narrower the 
bandwidth is, the smaller is the propensity of the excess 
carrier to delocalize. 

A common feature evident in results for charged chains 
in vacuum is appreciable overall energy shifts with re- 
spect to neutral chain MO levels: to higher energies for 
excess electrons and to lower energies for excess holes. 



These shifts reflect simple electrostatic effects (the elec- 
trostatic potential with respect to infinity) for chains of 
finite lengths. Understandably, shifts of such an origin 
should become much smaller for chains surrounded by 
media with high dielectric constants, which is indeed the 
case for our chains in the solvent environment (see also 
below in Sec. IIII C|) . A much more important and rele- 
vant difference in the MO structure of charged chains in 
vacuum and in solvent is instead in that solvated chains 
feature local electronic levels clearly separated from other 
states. In the case of negatively charged chains, Fig. [51 
it is the highest occupied molecular orbital (HOMO) ac- 
commodating an excess electron. In the case of posi- 
tively charged chains, Fig. it is the lowest unoccu- 
pied molecular orbital (LUMO) accommodating an ex- 
cess hole. Comparison of panels (d) to (b) and, espe- 
cially, panels (h) to (f) in these figures clearly demon- 
strates the formation of local levels - as would be caused 
by the appropriate polarization of the environment. The 
appearance of localized electronic states in self-consistent 
"potential wells" is one of the hallmarks of the strong po- 
laronic effect [ol [Tlj|. 

The number of localized electronic states appearing 
due to the polaronic effect should depend on the nature of 
the system and the effective potential well. For electron- 
lattice polarons in CPs, local levels are known [l|, Q to 
split off from both conduction and valence bands. The 
short-range character of the effective potential well in 
that case, however, limits the number of such states (in 
traditional two-band Peierls dielectric models |f|, there 
would be only two local levels: one close to the conduc- 
tion and one to the valence band). In the case of the 
long-range interaction with the polarizable environment, 
the effective potential well has the long-range Coulomb 
behavior expected to result, for an infinite system, in a 
multitude of localized states of various spatial symme- 
tries converging to the onset of the delocalized contin- 
uum. This is a standard picture for traditional "single- 
band" 3D polarons @ , and we have specifically discussed 
its realization for ID systems immersed in the 3D polar 
medium 0,0|- 

Restricting a more detailed illustration now to the neg- 
atively charged CiooH^" chains, Fig. [8] compares the spa- 
tial structure of HOMOs and several other MOs (from a 
list ordered in orbital energies) in vacuum and in the sol- 
vent environment as retrieved from the Gauss View 4 (4lj 
rendition of Gaussian 03 outputs. Here we use results 
only for RG2 geometry displaying better distinguishable 
localized states. (The figure of course also illustrates the 
lifting of spin and spatial-orientation degeneracies.) As 
expected, a stark difference is seen between the HOMOs 
in vacuum and solvent cases: the former shows a state de- 
localized over the chain length, while the latter indeed ex- 
hibits a well-localized character. In addition to this, how- 
ever, we also find in the right column of Fig. [8] that the 
solvation-induced negative polaron features various occu- 
pied localized states in the valence band region of energies 
(nothing like this happens for valence band states in the 
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FIG. 8: Spatial distribution of several MOs for CiooH^" chains with RG2 geometry in vacuum (left column) and in solvent 
(right column). The spin state of orbitals is indicated by arrows, while the spatial orientation by letters "x" and "y" (chains 
are oriented along the z-axis). Orbital energies in eV are shown to the right of the columns. Compared MOs around HOMOs 
are sequentially selected from a list ordered in orbital energies. 



vacuum case). This observation confirms a reorganiza- 
tion of valence band states apparent in the corresponding 
electronic structure of FigJBfh). Thus the polaron reveals 
a many-electron structure, more complex than invoked in 
simplified single-particle models [!, 0, H, H, 0] ■ 

Even in RG2 geometry it is difficult though to visu- 
ally identify the LUMO in the right column of Fig. [5] 
as an unoccupied localized state expected to be formed 
in the effective polarization potential well 0, 0] ■ While 
that state wave function evidently exhibits the right odd 
symmetry with respect to the inversion of the chain axis, 
the length of the chain appears insufficient to establish a 
localized character of the state with certainty. Likewise, 
within the length of our oligomer we cannot establish the 
onset of what would become a continuum of truly delo- 
calized states in the limit of infinite chains. 

Similar considerations of details of the spatial behavior 
of MOs can be done for positively charged chains which 
would then correlate with signatures of the electronic 
structure in Fig. [7J 



C. Energetics 




20 40 60 80 100 20 40 60 80 100 
N 



FIG. 9: The negative of electron affinities, — EA, (square 
symbols) and energies -Elumo of neutral chains (triangles) 
as functions of the number of carbons in polyyne oligomers 
CatH2 with rigid geometries. Vacuum data points are con- 
nected by dashed lines, data for chains in solvent by solid 
lines. Panel (a) displays results for chains with RG1 geome- 
try, panel (b) corresponds to chains with RG2 geometry. 



An important energetic quantity in discussions of the 
polaronic effect is the polaron binding energy referring 
essentially to the gain in the total energy of the system 
taking place upon self-localization of the excess charge 
carrier [l], 0, Calculations of the binding energy 

of solvation-induced polarons on ID structures are very 
straightforward in simplistic models [!, 0, H, [E 0] ■ Un- 
fortunately, we do not feel we have a reliable straightfor- 
ward procedure for determination of the binding energy 
of solvation-induced polarons in our current ab initio cal- 
culations. One of the reasons for this is that we cannot 
control the state of the polarization of the solvent as it 
would be easily done by prescribing desired patterns of 



atomic displacements in the case of electron-lattice po- 
larons. Some considerations can, however, be put for- 
ward based on direct outputs of Gaussian 03 computa- 
tions. 

Figure [9] illustrates such considerations for excess elec- 
trons. One very well-defined and physically meaningful 
quantity characterizing accommodation of an excess elec- 
tron in the computational outputs is the electron affinity 



corresponding to the difference of total system energies 
for neutral (superscript "0" ) and negatively charged ( "-" ) 



8 



chains. Another readily available energy is the energy 
^lumo °f the LUMO in the neutral system. These en- 
ergies are compared in Fig. [9] as — EA vs -Elumo f° r var " 
ious conditions as functions of the number N of carbons 
in Cj\rH2. The rationale for this comparison is that in 
the limit of infinite chains, N — > op, one might hope, in 
the spirit of Koopman's theorem [28|, [4^], that E® VMO 
would represent the extra system energy upon an addi- 
tion of a single electron into the lowest-energy delocalized 
state with vanishing excess charge density. The negative 
of the electron affinity, — EA, on the other hand, repre- 
sents the extra system energy upon addition of a single 
electron when the system is allowed to completely reorga- 
nize itself to achieve the minimum possible total energy. 
Both compared "extra energies" are of course negative 
in our data. We remind the reader that in this section 
we discuss system with rigidly positioned carbons so that 
reorganization here is limited only to the states of other 
electrons and the polarization of the environment. 

A striking difference is observed in Fig. [5] when re- 
sults for chains in vacuum are compared to chains in 
solvent. In the former case the data corresponds to 



-EA > E\ 



LUMO' 



-Extra [< > 



while in the latter the order of these 
energies is reversed: — EA < E® UMO , clearly showing 
the significance of screening and reorganization effects 
taking place in the solvent environment. Vacuum data 
points are evidently very far from displaying a conver- 
gence to the infinite chain limit, and a larger steepness 
of the vacuum EA as a function of N is indicative of 
the magnitude of the Coulomb effects involved for fi- 
nite excess charge densities in oligomers. While data 
points for chains in solvent, of course, also exhibit an 
iV-dependence, this dependence is shallower than in the 
vacuum case, and for RG2 geometry (panel (b)) displays 
a nearly converged behavior for largest N studied. As we 
discussed in Sec lIII Al a self-consistent polaron formation 
is achieved in RG2 within our range of chain lengths. 
Considering the positive energetic "improvement" 
EA in the solvent environment, its N = 100 
magnitudes in Fig. [9] are equal to approximately 0.1 and 
0.11 eV for RG1 and RG2 geometries, respectively. If one 
were to take this energetic improvement in the infinite- 
TV limit as corresponding to the "true" solvation-induced 
polaron binding energy, then the obtained numbers ~ 0.1 
eV would provide estimates of the polaron binding. 

An analogous consideration for excess holes involves a 
comparison of the ionization potential 

IP = E+ t - £° ot 

and the negative — -Ejjomo °f ^ ne energy of the HOMO in 
the neutral chain. This comparison of now positive extra 
system energies is illustrated in Fig. |TD] The figure dis- 
plays qualitative patterns similar to observed above for 
excess electrons (while, of course, also showing signs of 
CCS-breaking). Particularly, one notes while for chains 
in vacuum one has IP > — -Ehomc ^ ne relationship re- 
verses for chains in solvent: IP < — -Ehomo- The N = 80 
values of the energetic "improvement" — -Ejjomo — ^ 



7.4- 
7.2. 

^7.0 

> 



^6.8 
>> 
M 

6.6 

a 

W 

6.4 
6.2 
6.0 



(a) 




(b) 



9.2 
9.0 
8.8 
8.6 
8.4 
8.2 
8.0 
7.8 



20 



40 



60 



80 20 
N 



40 



60 



80 



IP, (square symbols) and 



J?homoi °f neutral chains 



FIG. 10: Ionization potentials 
the negative of HOMO energies, 
(triangles) as functions of the number of carbons in CjvH2 
with rigid geometries. Vacuum data points are connected by 
dashed lines, data for chains in solvent by solid lines. Panel 
(a) displays results for chains with RG1 geometry, panel (b) 
corresponds to chains with RG2 geometry. 



in the solvent environment are approximately 0.1 and 
0.08 eV for RG1 and RG2 geometries, respectively, which 
could again suggest polaron bindings close to 0.1 eV. 

It appears, however, strange that these naive estimates 
lead to relatively small differences (and even to a rever- 
sal) in polaron bindings for RG1 and RG2 geometries - 
while the spatial localization has been found much more 
significant in RG2, Sec. IIII Al We believe this indicates 
that the actual polaron binding energies would be under- 
estimated in discussed energetic comparisons. It is quite 
likely that El VMO and — -^homo °-0 11011 represent fully 
the lowest extra energy of systems where an excess de- 
localized electron or an excess delocalized hole would in- 
teract with a non-uniform, charge density wave (CDW), 
background of the polymer chain. 

Bond-centric CDWs are known to form on neutral 
chains upon their dimerization with more electronic den- 
sity corresponding to shorter (C = C) and less electronic 
density to longer (C — C) bonds. Such CDWs have been 
widely discussed in the context of CPs as Peierls insula- 
tors [Ij . Being a non- vanishing distribution of the electric 
charge, a CDW is expected to polarize the solvent, which 
is indeed how we interpret the data shown in Fig. 111! The 
figure displays the electrostatic energy U po i stored in the 
solvent polarization as retrieved from outputs of Gaus- 
sian 03 computations. It is evident that this energy does 
not vanish for neutral chains and in fact grows linearly 
with the number of carbons N, as one would expect for a 
CDW. Since the amplitude of CDW is enhanced in neu- 
tral chains with RG2 geometry, they do exhibit a steeper 
growth with N. 

Figure [IT] also shows the evolution of U po \ for neg- 
atively charged CatH^~ chains. Especially illuminating 
here are the results for chains with RG2 geometry ex- 
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FIG. 11: The electrostatic energy stored in the solvent polar- 
ization as a function of the number N of carbons in CjvEb. 
Triangular data points are for neutral CnYQ, squares for 
charged CjvH^" chains. Dashed lines connect the data for 
RG1 chain geometry, solid lines for RG2. 



hibiting a non-monotonic behavior of U po \ as a function 
of N. E/poi first decreases with N at smaller N reflecting 
the decrease of the electric field of the excess charge upon 
its spreading over longer oligomers. U po \ then reaches its 
minimum in the region of N corresponding to the for- 
mation of the self-consistent localized polaron pattern. 
With a further increase of N, U po \ starts growing again 
with the same slope that is exhibited by neutral chains 
reflecting the interaction of the CDW with the solvent. 
The converged large- N difference in U po \ of charged and 
neutral chains with RG2 geometry is evidently quite sub- 
stantial, close to 0.9 eV. While the minimum of U po \ is not 
achieved for negative chains with RG1 geometry within 
the range of chain lengths studied, one could guess from 
the data of Fig. [11] that the large- N difference in U po \ of 
charged and neutral chains in that case would probably 
be roughly half of the value for RG2 geometry. This ratio 
would look like more in line with charge localization pat- 
terns observed for RG1 and RG2. The very magnitude 
of the "excess" polarization energy just discussed may 
be indicative of the polaron binding energies apprecia- 
bly larger than 0.1 eV. This speculative conclusion could 
be verified more directly if we were able to perform ab 
initio computations for charged chains with the polariza- 
tion state of the solvent corresponding to that for neutral 
chains. 



IV. SOLVATION AND BOND ALTERNATION 
PATTERNS 

We now turn our attention to the results obtained with 
full geometrical optimizations of underlying carbon lat- 
tices. A very convenient and well-known [H [2a. l27j way to 
characterize dimerized polymeric structures is via bond- 



length alternation (BLA) patterns 

*»=(-l) W (In-ln+l), 

where l n is the length of the nth carbon-carbon bond. 
In the infinite dimerized neutral structure, this pattern 
would be uniform, that is, independent of the spatial 
bond position n. 

Figure [T^] compares BLAs computed for neutral and 
negatively charged oligomers C40H2, C80H2 and C100H2. 
This comparison allows one to visually evaluate the mag- 
nitude of oligomer end effects. It is clear that a uniform, 
infinite-polymer, behavior is practically achieved over a 
long central-part stretch of the neutral CiooH!]. There 
is very little difference in BLAs of neutral chains in vac- 
uum and in solvent; a small increase of 8 n in the solvent 
is consistent with our expectation for a CDW interacting 
with the environment. 

Much more significant variations are observed for 
charged oligomers. As was discussed at length for 
electron-lattice polarons in CPs [l|, l2| , the magnitude of 
BLAs would exhibit a dip in the spatial region contain- 
ing an excess particle, where the lattice becomes "less 
dimerized" . This expected decrease of 6 n is evident in all 
panels of Fig. [T2j The character of its evolution with the 
oligomer length, however, is quite different for chains in 
vacuum and in solvent. The vacuum data just shows a 
progressive reduction of the dip in S n for longer oligomers. 
In the absence of data for much longer chains, all we can 
conclude is that a very flat pattern of 5 n in the middle 
of CiooH^ does not suggest yet an incipient formation of 
a self- localized pattern characteristic of polarons. This 
is reminiscent of the apparent failure of DFT computa- 
tions to detect electron-lattice polarons in computation- 
ally tractable oligomers discussed in Sec. [J- with our re- 
sults now extended to hybrid B3LYP computations. On 
the contrary, the data derived for chains in solvent quite 
clearly indicates a nearly self-consistent localized pattern 
of S n on CiooH^, although, of course, longer oligomers 
would be needed to reveal a truly self-consistent polaron 
size. 

These observations correlate well with the results ob- 
tained for excess charge densities on CiooH^ that are 
displayed in panel (a) of Fig. [L3j The panel compares 
the excess density on chains with the uniform BLA of 
neutral chains (RG1 geometry) and on chains where the 
relaxation of the BLA due to the excess charge has taken 
place. One can see that while the relaxation of the BLAs 
does lead to spatial contractions of charge densities both 
in vacuum and in solvent, the delocalized (over the whole 
oligomer) nature of the excess charge remains intact in 
the vacuum case. Excess charge localization due to solva- 
tion, on the other hand, gets evidently enhanced by the 
localized lattice response seen in Fig. H2T c). One could 
probably say that solvation and lattice relaxation here 
act synergistically to affect excess charge self-trapping. 
With our B3LYP computations, the effect of the solva- 
tion on the excess charge distribution is evidently much 
more pronounced. A complementary information is pro- 
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FIG. 12: Carbon-carbon bond-length alternation patterns for (a) C40H2 (b) C80H2 and (c) C100H2. Patterns for neutral chains 
are indicated with symbol (0), patterns for negatively charged chains with (-1). Dashed lines show results for chains in vacuum, 
solid lines for chains in the solvent environment. 
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FIG. 13: Comparison of the solvation and BLA effects for 
C100H2. (a) Excess charge distributions in vacuum and in the 
solvent for the fully optimized lattice geometry (solid lines) 
and for the RG1 geometry (dashed lines), (b) The corre- 
sponding MO energy levels in the solvent for: (bl) Optimized 
neutral chain, (b2) Charged chain in RG1 geometry, (b3) Op- 
timized charged chain. The arrows show HOMO levels. 



vided by panel (b) of Fig. [13] showing the effects of sol- 
vation and BLA relaxation on the electronic structure of 
MO levels. Particularly relevant here are the intragap 
levels. Interestingly, for the energies of these levels ef- 
fects due to each of the mechanisms appear similar in 
magnitude. This suggests that one has to be careful in 
evaluation of self-localization per se based just on the 
appearance of the MO levels in their overall structure. 
The spatial behavior of individual MOs as well as the 
total spatial distribution of the excess charge need to be 
properly analyzed. 



A. Hartree-Fock results 

As mentioned in the Introduction, in this paper we 
do not pursue an extensive comparison of results ob- 
tained with various ab initio methods. To emphasize the 
level-of-theory dependence of at least the quantitative as- 
pects, we, however, find it pertinent to provide a glimpse 
at the nature of results derived with Hartee-Fock (HF) 
computations. HF calculations for conjugated polymeric 
structures have been discussed [HI, |27| | to result both in 
overestimated energy gaps and magnitudes of BLAs. At 
the same time, unlike pure-DFT computations, they have 
also been shown to yield self-localized polaronic states for 
excess charge carriers [22|, [23|, [24], HH • Our findings for 
polyynic chains agree well with these known trends. 

Figure [Til illustrates our restricted-HF results for neu- 
tral and negatively charged C40H2 chains both in vac- 
uum and in the solvent environment; they are consistent 
with a limited set of data we have also computed for 
longer chains. When compared to B3LYP-derived re- 
sults discussed earlier in this paper and in Ref. [l5j|, one 
immediately notices an appreciably higher degree of ex- 
cess charge density localization and larger magnitudes of 
BLAs obtained with the HF. Reassuringly, the qualita- 
tive effect of the solvent environment remains intact: it 
does enhance the excess charge localization (panel (a) of 
Fig. [T4|) and leads to a shorter and more pronounced asso- 
ciated dip in the optimized BLA pattern (panel (b)). Dif- 
ferently from the B3LYP-data, however, solvation here 
does not appear to act as a "primary" source. 

In agreement with previous HF studies, we do find the 
excess charge self-localized already on chains in vacuum. 
One might be tempted to say - as expected for tradi- 
tional electron-lattice polarons - that the calculated self- 
trapping in this case is caused by the interaction with 
displacements of carbon atoms manifested in the relaxed 
BLA pattern in Fig. [Tjjb). The comparison made in 
Fig. [LlTa). however, does not seem to support such a 
viewpoint. Specifically, we compare in that panel the 
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FIG. 14: Results of restricted-HF computations for C40H2 
oligomers, (a) Excess charge density distribution on C4oH^" 
in vacuum and in solvent. Dashed lines show the distribu- 
tion on chains with a uniform BLA pattern of the neutral 
system, solid lines are for chains with relaxed (optimal) ge- 
ometries, (b) Optimal BLA patterns for neutral (0) and neg- 
atively charged (-1) oligomers. Dashed lines are for chains in 
vacuum, solid for chains in solvent. 



excess charge densities obtained for fully relaxed geome- 
tries and for a rigid geometry featuring a uniform BLA 
pattern. The latter corresponds to the BLA derived for 
neutral chains in their middle portion, Fig. RHT b). De- 
spite the uniform lattice background, the excess charge 
localization evidently takes place, and with minimal vari- 
ations from the fully geometrically relaxed chains. More- 
over, the lattice relaxation, while indeed decreasing the 
total system energy, counter-intuitively results here in 
slightly smaller localization. (To compare, in B3LYP re- 
sults geometrical relaxation would always lead to larger 
spatial charge localization.) Thus, it is rather the treat- 
ment of electron-electron interactions within the frame- 
work of HF computations, that appears to effect the ex- 
cess charge localization observed in our HF data. In this 
paper, we forgo a discussion of the physical validity of 
these results that would require a more focused study. 



V. DISCUSSION 

Self-localization of excess charge carriers (electrons or 
holes) into polaronic states on ID semiconductor struc- 
tures is an interesting effect that may be consequential 
for transport, optical, and electrochemical properties of 
these systems. The basic nature of the self-trapping 
is easily revealed in a single-particle picture within a 
standard ID continuum adiabatic framework, where a 
localized carrier wave-function tp(x) is found as corre- 
sponding to the ground-state solution of the non-linear 
Schrodinger equation: 



h 2 d 2 jj(x) 

2m dx 2 



dx'V(x - x) \*p(x')\ 2 tp(x) = Eifj(x) 



Here x is the coordinate along the structure aixs, m the 
effective band mass of the carrier, and V(x) the effective 
self-interaction mediated by another subsystem. In the 
case of a short-range electron-phonon mediation, the self- 
interaction may be approximated as local: V(x) — gS(x) 
leading then to the exact solution ip{x) cx sech(gmx/2h 2 ) 
for a continuum ID polaron. This is the result widely 
known after the pioneering contributions of Rashba [43[ 
and Holstein [2(|. Another type of mediation arises due 
to the long-range polarization interaction of an excess 
charge carrier with a surrounding polar medium. In that 
situation, the effective V(x) behaves as l/\x\ at large dis- 
tances, while the short-range behavior depends on spe- 
cific system details such as, e.g., the actual transverse 
charge density distribution and the geometry of the di- 
electric screening [1, [3, @] • Numerical studies of this case 
d, [|| indicate that the binding energy of the resulting po- 
larons could reach an appreciable fraction of the binding 
energy of the corresponding primary optical excitations, 
Wannier-Mott excitons. 

While illuminating the basic physics of the self- 
trapping, simplified single-particle descriptions have the 
drawback of not explicitly including valence-band elec- 
trons, whose interactions and reorganization may sub- 
stantially affect the outcomes. Numerous ab initio calcu- 
lations have been dedicated to effects of electron-electron 
interactions on electron-lattice polarons in CPs . In this 
pap er we have extended our first model ab initio study 
[lj| intended to address these effects in the formation 
of ID polarons due to the interaction with a surround- 
ing polar solvent. As our study is performed within the 
framework of the polarizablc continuum model (PCM), 
it could also be thought of as a new avenue for applica- 
tions of quantum-mechanical continuum solvation models 

MM- 

By studying and comparing various model systems 
based on long CatH2 oligomers with the polyynic struc- 
ture, we have been able to successfully demonstrate 
the self-consistent formation of ID large electron- and 
hole-polarons entirely due to solvation. Reorganization 
of valence-band electrons and a more complex many- 
electron structure of the resulting polarons have been 
also explicitly illustrated. Full geometrical optimization 
of the underlying carbon lattice has shown that localized 
bond-length alternation patterns may act synergistically 
with the solvent reorganization to even further increase 
the degree of excess charge self-localization. 

It should be reiterated that, while providing clear con- 
ceptual demonstrations and a seemingly appealing physi- 
cal picture of solvation-induced self-localization, our cur- 
rent study is limited both in terms of its scope and direct 
quantitative applicability to real systems. Among other 
things, we have already discussed and illustrated a rela- 
tively strong dependence of the results on the level of ab 
initio theory used, a common feature of such computa- 
tions [2^, [27| . DFT computations with hybrid exchange- 
correlation functionals may be "appropriate" engines for 
the problem at hand but a more comprehensive compar- 
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ative study is required before conclusions can be drawn. 
The PCM framework as used in this paper may overesti- 
mate the effects of solvation on self-localization as it does 
not take a proper account of the frequency dependence of 
the actual solvent dielectric functions 

HIE IE El- This 

issue could possibly be analyzed within newly developed 
time-dependent PCM schemes 0, 5E1- As the compu- 
tational demand on calculations of ID large polarons, 
especially in application to more complex semiconductor 
structures, may turn out to be excessive, one should also 
be looking for more efficient computational approaches. 



It would be interesting to see if the scaling description 
of solvent effects (46| may be applied to the problem of 
charge-carrier self-localization on ID semiconductors. 
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